Thermodynamic anomalies in a lattice model of water: Solvation properties 
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We investigate a lattice-fluid model of water, defined on a 3-dimensional body-centered cubic 
lattice. Model molecules possess a tetrahedral symmetry, with four equivalent bonding arms. The 
model is similar to the one proposed by Roberts and Debenedetti [J. Chem. Phys. 105, 658 
(1996)], simplified by removing distinction between "donors" and "acceptors". We focus on solvation 
properties, mainly as far as an ideally inert (hydrophobic) solute is concerned. As in our previous 
analysis, devoted to neat water [J. Chem. Phys. 121, 11856 (2004)], we make use of a generalized 
first order approximation on a tetrahedral cluster. We show that the model exhibits quite a coherent 
picture of water thermodynamics, reproducing qualitatively several anomalous properties observed 
both in pure water and in solutions of hydrophobic solutes. As far as supercooled liquid water is 
concerned, the model is consistent with the second critical point scenario. 
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I. INTRODUCTION 

From the experimental point of view, water is known to 
exhibit several thermodynamic anomaliesii2i£. Contrary 
to most fluids, at ordinary pressure the solid phase (ice) 
is less dense than the liquid phase. The latter displays 
a temperature of maximum density at constant pressure, 
while both isothermal compressibility and isobaric heat 
capacity display a minimum as a function of temperature. 
Moreover, heat capacity is on average much larger than 
usual. Anomalous properties of neat water have been 
studied for long, but much interest has been devoted also 
to unusual properties of water as a solvent, in particular 
for nonpolar (hydrophobic) chemical species^* 5 *^. Inser- 
tion of a nonpolar solute molecule in water is character- 
ized by positive solvation Gibbs free energy (it is ther- 
modynamically unfavored), negative solvation enthalpy 
(it is energetically favored), negative solvation entropy 
(it has an ordering effect), and large positive solvation 
heat capacitji&S. More precisely, for prototype hydropho- 
bic species (that is, for instance, noble gases), solvation 
entropies and enthalpies, which are negative at room 
temperature, increase upon increasing temperature, and 
eventually become positive. These properties define the 
so-called hydrophobic effect, whose importance in biolog- 
ical processes, such as protein folding, has been empha- 
sized in the latest yearsiS. 

From the theoretical point of view, one can relate the 
anomalous properties of neat water to the formation of 
a large amount of hydrogen bonds, because of peculiar 
features of water moleculesii*i£. The same physics is be- 
lieved to underly the hydrophobic effec t 4 ! 5 ' 13 , but a com- 
prehensive theory which explains all of these phenomena 
has not been developed yet. A possible way of inves- 
tigation consists of computer simulation o 14 : 15 : 16 , based 
on very detailed (but still phcnomenological) interaction 
potentials. In this way, quite a high level of accuracy 
in describing water thermodynamics has been achieved. 
Nevertheless, simulations are generally limited by the 



large computational effort required, while microscopic 
physical mechanisms are sometimes hidden by a large 
number of model details. A complementary approach 
involves investigations of simplified modelsSiiiiSii^. Al- 
though quantitative accuracy is sometimes poor, this ap- 
proach usually allows more detailed analysis, in a wide 
range of thermodynamic conditions, with relatively low 
computational cost. One of these attempts is based on 
the application of scaled-particle theorjiSS to hydropho- 
bic hydration^. A recent descendant of scaled-particle 
theory is the information theory approach by Pratt and 
coworkers^, based on previous knowledge of water prop- 
erties, such as the pair correlation function, which can 
be obtained by experiments or by simulations. Accord- 
ing to the cited studies, the hydrophobic effect would 
result mostly from small size of water molecules, and not 
from water structuring by the solute, as in the classical 
view4. Such effect, though existing, would be scarcely 
relevant for a description of the hydrophobic effect. The 
simplified molecular thermodynamic theory of Ref. ^3 
is basically in agreement with this conclusion. Differ- 
ent theories stress that the large positive heat capacity 
variation, observed upon insertion of apolar solutes into 
water, can only arise from cooperativity, that is, from 
induced ordering of water molecules, so that a theory 
of the hydrophobic effect should be based on a descrip- 
tion of this phenomenon. This position is supported by 
the results of the 2-dimensional "Mercedes Benz" model, 
first introduced by Ben-Nairn in 1971»2i, and extensively 
investigated by Dill and coworkers in the latest yearst 
Contrary to the previously mentioned approaches, the 
Mercedes Benz model, though involving high simplifica- 
tions, is based on well defined microscopic interactions, 
that is, on an energy function, without previous knowl- 
edge of water properties. One important reason to do so 
is the need of modelling water in a computationally con- 
venient way, in order to investigations on complex sys- 
tems such as biomolecules, for which water plays a key 
role. An even more simplified approachi&S 4 ^, which 
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nevertheless can in principle satisfy this criterion, re- 
lies on the long standing tradition of lattice fluid mod- 
els. As far as neat water is concerned, several differ- 
ent models, both in 2 26 i 27 i 28 i 29 i 30 i 31 i 32 i 33 and 3 dimen- 
siona 34 i 35 i 36 i 37 i 38 i 39 i 40 i 41 have been investigated, some of 
which are variations of the early model proposed by Bell 
in 197224. One of them is the 3-dimensional model by 
Roberts and Debenedett i 40 i 42 i 43 , defined on the body- 
centered cubic lattice. In this model, water molecules 
possess four bonding arms (2 donors and 2 acceptors) ar- 
ranged in a tetrahedral symmetry. Hydrogen (H) bond 
formation requires that 2 nearest neighbor molecules 
point respectively a donor and an acceptor towards each 
other. A number of nonbonding configurations is al- 
lowed, to account for H bond directionality. Bond weak- 
ening occurs (the bond energy is increased of some frac- 
tion) whenever a third molecule is placed near a formed 
bond. The latter feature basically mimics the fact that 
too closely packed water molecules disfavor H bonding. 
Let us notice that, while bonding properties are equiva- 
lent to those of the Bell model 3 ^, the weakening criterion 
is different. The Roberts-Debenedetti model is quite ap- 
pealing in that it has been shown to predict not only 
some of real water thermodynamic anomalies, such as 
the temperature of maximum density, but also a liquid- 
liquid phase separation in the supercooled region, and 
a second critical point. Nevertheless, the distinction be- 
tween donors and acceptors is likely to be not so crucial to 
the physics of water. Therefore, in a previous paper™, 
we have investigated a simplified version of the model 
(without donor/acceptor distinction), showing that the 
same basic properties could be reproduced. Here we ex- 
tend the simplified model to deal with aqueous solutions, 
working out solvation thermodynamics for an inert (ap- 
olar) solute. Our purpose is to verify whether the model 
is able to reproduce also the main features of hydropho- 
bicity. This analysis might be interesting also in view of 
investigations on mixtures of water with more complex 
chemical species, such as polymers. We shall carry out 
the analysis by means of a generalized first-order approx- 
imation on a tetrahedral cluster, which has been verified 
to be quite accurate for the neat water model 4 ^. 



II. THE MODEL 

Let us introduce the model. Molecules are placed on 
the sites of a body-centered cubic lattice, whose structure 
is sketched in Fig. ^ A site may be empty or occupied 
by a water molecule (w) or by a solute molecule (s). An 
attractive potential energy — < is assigned to any 
pair of nearest neighbor (NN) sites occupied by molecules 
of species x, y, where x and y can take the values w, s. 
This is the ordinary Van der Waals contribution. Water 
molecules possess four equivalent arms that can form H 
bonds, arranged in a tetrahedral symmetry, so that they 
can point towards 4 out of 8 NNs of a given site. There 
is no distinction between donors and acceptors, so that 
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FIG. 1: Two conventional cells of the body-centered cubic 
lattice: A, B, C, D denote 4 interpenetrating face-centered cu- 
bic sublattices. 

a H bond is formed whenever two NN molecules have 
a bonding arm pointing to each other, yielding an en- 
ergy — rj < 0. It is easily seen that water molecules have 
only 2 different configurations in which they can form 
H bonds (see Fig. EJ) . We assume that w more configura- 
tions are allowed, in which water molecules cannot form 
bonds (w is related to the bond-breaking entropy). More- 
over, we assign an energy increase rjc x /6, with c x £ [0, 1], 
for each of the 6 sites closest to a formed bond (i.e., 3 out 
of 6 second neighbors of each bonded molecule), occupied 
by an x molecule. A bond surrounded by 6 molecules of 
species x contributes an energy — r;(l — c x ). As far as 
water molecule are concerned, the weakening parameter 
mainly accounts for the fact that H bonds are most fa- 
vorably formed when water molecules are located at a 
certain distance, larger than the optimal Van der Waals 
distance. Therefore, if too many molecules are present, 
the average distance among them is decreased, and hy- 
drogen bonds are (on average) weakened. Moreover, the 
presence of an external molecule may perturb the elec- 
tronic density, resulting in a lowered H bond strength as 
well. A weakening parameter for the solute (c s ) takes 
into account possible perturbation effects for a generic 
chemical species, even if in the following we shall mainly 
consider an ideally inert solute with c 5 = 0. 

The hamiltonian of the system can be written as a sum 
over irregular tetrahedra, whose vertices lie on 4 different 
face-centered cubic sublattices, shown in Fig. ^ O ne of 
such tetrahedra is shown in Fig. E^a). We have 

^ = g ^ ] Tl-iaipi^is > (1) 
(a,/3 )7 ,<5) 

where TLijki is a contribution which will be referred to as 
tetrahedron hamiltonian, and the subscripts i a ,i0,ij,is 
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FIG. 2: Two model molecules forming a H bond. The lower 
molecule is in the i = 1 configuration, the upper one is in the 
i = 2 configuration. 




FIG. 3: (a) Basic cluster (irregular tetrahedron): A,B,C,D 
denote sites in the 4 corresponding sublattices. AB, BC, CD, 
and DA are NN pairs; AC and BD are second neighbor pairs, 
(b) Husimi tree structure corresponding to the generalized 
first order approximation on the tetrahedron. 



label site configurations for the 4 vertices a, /3, 7, S, re- 
spectively. Possible site configurations are: "empty" 
(i = 0), "bonding water" (site occupied by a water 
molecule in one of the 2 orientations which can form 



TABLE I: Site configuration labels (i), with corresponding 
multiplicities (w;) and occupation numbers for water (n„ t i = 
1) and solute (n Sj i = 1). 
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bonds: i = 1,2; see Fig- EJ j "nonbonding water" (site 
occupied by a water molecule in one of the w orienta- 
tions which cannot form bonds: i = 3), "solute" (site 
occupied by a solute molecule: i — 4). Assuming that 
{j,k), (k,l), and refer to NN pair configura- 
tions, the tetrahedron hamiltonian reads 

'Hijki = Hijki + Hjku + Hkuj + Hujk, (2) 

where 

u h (^ n x,k+n x ,l\ ,„s 

Hijki = -^x y nx,in y ,j - yhij II- c x - I , (3) 

Ux,i is an occupation variable for the x species, defined 
as in Tab. while hij are bond variables, defined as 
hij — 1 if the pair configuration forms a H bond, 

and hij = otherwise. Here and in the following, we un- 
derstand that repeated x, y indices are to be summed over 
their possible values w, s. Let us also assume that i, j, k, I 
(in this order) denote configurations of sites placed on, 
say, A, B, C, D sublattices respectively. If A, B, C, D are 
defined as in Fig. ^ we can define hij = 1 if i = 1 and 
j = 2, and h^ = otherwise. With the above assump- 
tions, the tetrahedron hamiltonian is independent of the 
orientation, that is of the arrangement of A, B, C, D on 
its vertices. Let us notice that both Van der Waals 
(— e xy n Xi in y j) and H bond energies (— ry/iy), which are 
2-body terms, are split among 6 tetrahedra, whence the 
1/6 prefactor in Eq. Q). On the contrary, the 3-body 
weakening terms (rjc x hijn x .k/6) are split between 2 tetra- 
hedra, thus the 1/6 factor is absorbed in the prefactor, 
while a 1/2 factor is left in the tetrahedron hamiltonian. 

III. FIRST ORDER APPROXIMATION 

We perform the investigation by means of a generalized 
first order approximation on a tetrahedral cluster. In the 
previous paper—, we have introduced the approximation 
in the variational approach—, as a particular choice of the 
largest clusters left in the entropy expansion (basic clus- 
ters). Such a choice, sometimes denoted as cluster-site 
approximation^ (the only clusters to be taken into ac- 
count in the expansion are basic clusters and single sites), 
has not only the advantage of high simplicity, but also of 
relative accuracy, which has been recognized for different 
models 31 ' 45 . The basic clusters are a number of irregu- 
lar tetrahedra, namely, 4 out of 24 tetrahedra sharing a 
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given site, as sketched in Fig. [31 This choice turns out to 
coincide with the (generalized) first order approximation 
(on a tetrahedron), also equivalent to the exact calcula- 
tion on a Husimi latticed, whose (tetrahedral) building 
blocks are just arranged as in Fig. OJb). In the present 
paper, we employ the latter approach, which, not hav- 
ing to treat cluster probability distributions explicitly, is 
numerically more convenient. 

Husimi lattice thermodynamics can be studied exactly, 
in a numerical way, by solving a suitable recursion rela- 
tion 46 , since the system is locally treelike. First of all, 
we have to choose, as a Husimi lattice hamiltonian, the 
following expression 

W = ^2 T~Li a ipiyi6i ( 4 ) 

obtained from Eq. JIJ by understanding that the sum 
runs over tetrahedra in the treelike system only, and re- 
moving the 1/6 prefactor. Let us notice that the latter 
change is required, in order to obtain equal internal en- 
ergy densities (internal energies per site) for the Husimi 
lattice and the ordinary lattice system. If we denote the 
tetrahedron configuration probability by Pijki, and as- 
sume that it is equal on every tetrahedron, the internal 
energy density can be written as 

4 4 4 4 

U = X! Wi X! U 'J X! Wk X! w lPijkl7~Li]kU (5) 
j=0 j=0 k=0 1=0 

where Wi is the multiplicity of the i-th site configuration, 
equal to w for the nonbonding configuration of water 
(i = 3) and equal to 1 otherwise (see Tab-QJ. One then 
has to define partial partition functions (PPFs) in the fol- 
lowing way. Let us consider a single branch of a Husimi 
tree, made up of tetrahedral blocks, and a corresponding 
partial hamiltonian, obtained by Eq. with the sum 
restricted to tetrahedra in the branch. The PPF Qi is 
defined as a sum of the Boltzmann weights of the partial 
hamiltonian over the configurations of the branch minus 
the base site (that is why the PPF depends on a site con- 
figuration variable i). Working in the grand-canonical en- 
semble, as we are interested in, we also have to take into 
account chemical potential contributions, which is done 
by replacing the tetrahedron hamiltonian TLijki with 

~% n, n x .i + n Xl j + n x .k + n x .i 

rtijkl = rLijkl — Mx t j (o) 

where the usual convention on repeated indices holds. 
Let us notice that of course the PPF tends to infinity 
in the thermodynamic limit, that is for an "infinite gen- 
eration branch" , therefore it is convenient to define a 
normalized PPF qi (x Qi, for instance in such a way that 

4 

^2w t q l = l. (7) 

i=0 



In this way, qi represents the i-th configuration probabil- 
ity that the base site would have if it were not attached 
to any other branch. 

Let us now consider again the single branch of our 
Husimi tree. In the infinite generation limit, and in the 
hypothesis of a homogeneous system, the subbranches 
attached to the first tetrahedral block should be equiva- 
lent to the main one, so that one can write the recursion 
relation 

4 4 4 

qi = 2T 1 w i Y, Wk Y, u>ie-** w/r {WkQi? , (8) 
j=o k=a i=o 

where the sums run over configuration variables in the 
tetrahedron except i, T is the temperature expressed in 
energy units (whence entropy will be expressed in natu- 
ral units), and y is a normalization constant, imposed by 
Eq. J7J. The recursion relation can be iterated numeri- 
cally to determine a fixed point, representing the PPF of 
a branch whose base site lies in the bulk of the Husimi 
tree (generally denoted as Husimi lattice). Husimi lat- 
tice properties are equivalent to those obtained by the 
cluster-site approximation 46 . We can compute the site 
probability distribution p$, by considering the operation 
of attaching 4 equivalent branches to the given site. We 
obtain 

p t = z~\t, (9) 

where 

4 

z = Y,mqf (10) 

i=0 

provides normalization. We can also compute the tetra- 
hedron probability distribution, by considering the oper- 
ation of attaching 3 equivalent branches to each site of a 
given tetrahedron, yielding 

Pijkl oc e ~ n * jkl/T (qiqjqkqi f , (11) 

where of course the proportionality constant is deter- 
mined by normalization. From the knowledge of the 
tetrahedron probability distribution {pijki} one can com- 
pute the thermal average of every observable in the 
first order approximation, the internal energy density by 
Eq. Q and the grand-canonical free ener gy b y Eq. (10) 
in Ref. 41. According to Eq. (31) in Ref. U^L the latter 
can be also related to normalization constants as 

uj = -TQny - 21nz) , (12) 

where y is the normalization constant of the recursion re- 
lation (JSJ) and z is given by Eq. (|10|) . Finally, the entropy 
density can be computed as 



where u — [i x p x has formally the same expression as u in 
Eq. JSJ, with the tetrahedron hamiltonian TiijM replaced 
by Hijki ■ 
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FIG. 4: Pressure (P/rj) vs temperature (T/rf) phase diagram 
for pure water with e^/rj = 0.3, c w = 0.5, and w = 20. 
Thick solid lines denote (first order) phase transitions; circles 
are critical points. Thin (solid and dash-dotted) lines denote 
the spinodals and the Kauzmann line, respectively. Dotted 
straight lines denote constant temperature (T/77 = 0.34) and 
constant pressure {P/r) = 0.015) "cuts", corresponding to the 
solution phase diagrams, reported in Fig. |S| 



IV. RESULTS 

In order to investigate the model properties, we fix a 
set of parameters, as a result of several attempts. First 
of all, we take e ww /?/ = 0.3. This value of water- water 
Van der Waals interaction is equal to the one employed 
for the very detailed analysis by Roberts et al— , though 
a bit larger than previously employed by usai. Anyway, 
this choice accounts for the greater binding energy of hy- 
drogen bonds with respect to Van der Waals interactions. 
As far as the multiplicity of nonbonding water configu- 
rations is concerned, we set w = 20 (as in the previous 
work), which mimics the high directionality of hydrogen 
bonds. For neat water, it is necessary to set this parame- 
ter large enough to let anomalous properties appear, but 
further increase does not change qualitatively the phase 
diagram and the thermodynamic properties. As far as 
the weakening parameters are concerned, for water we 
choose c w = 0.5, which, given the other parameter val- 
ues, corresponds to a situation without a reentrant spin- 
odal, in agreement with most recent molecular dynamics 
results 3 ". 



A. Pure water properties 

Since the parameter choice is slightly different from the 
former work^i, we first reconsider neat water properties, 
taking the limit /i s — > — oo. This analysis is mainly meant 
to show that the model behavior still remains consistent 
with a "realistic" phase diagram, that is, with several ex- 
perimental evidences as well as simulation predictions. In 
fact, we have observed that relatively small variations in 
parameter values may give rise to quite dramatic changes 



in the phase behavior—, as observed also in other models 
of water—. The temperature-pressure phase diagram is 
reported in Fig. 21 Imposing homogeneity, our analysis 
includes both thermodynamically stable and metastable 
(supercooled) phases, even if stability is not investigated. 
Let us notice that pressure can be determined as P = —w, 
where u> is the grand-canonical potential per site, intro- 
duced in the previous section. We have assumed the vol- 
ume per site equal to unit, so that pressure is expressed 
in energy units. The appropriate order parameter is the 
density, i.e., the probability p w that a site is occupied 
by a water molecule, which can be evaluated from the 
formula 

4 

Px = »•,/),//....,. (14) 

i=0 

We find two different first order transition lines, ter- 
minating in two different critical points. The positively 
sloped line, at lower pressures, corresponds to the coex- 
istence of a very low density phase and a high density 
phase, and represents the ordinary vapor-liquid transi- 
tion. The other one, negatively sloped and placed at 
higher pressures, corresponds to coexistence between the 
high density liquid and and a lower density one. The 
related critical point may reasonably represent the so- 
called second critical point, which has been conjectured 
and observed in simulations, and of which also some ex- 
perimental evidences have been given. As one could ex- 
pect, the low density liquid turns out to be more hydro- 
gen bonded than the high density one. We find a density 
maximum as a function of temperature for liquid coex- 
isting with vapor (and at constant pressure as well), and 
the temperature of maximum density slightly decreases 
upon increasing pressure, as observed in experiments. 

We also report the liquid phase spinodals and the 
Kauzmann line. Details about the (semi-analytical) cal- 
culation of spinodals, which allows to determine density 
response functions as well, are given in the previous pa- 
per—. The limit of stability of the liquid phase (spin- 
odal) is the locus at which the metastable liquid ceases 
to be a minimum of (a variational form of) the free en- 
ergy, and becomes a saddle point. On the contrary, the 
Kauzmann line is the locus at which the liquid phase en- 
tropy vanishes, and corresponds to the ideal glass transi- 
tion. It can be easily determined numerically. As previ- 
ously mentioned, we do not observe a reentrance of the 
liquid-vapor spinodal in the positive pressure half-plane, 
which is actually a possibility of our model, for a dif- 
ferent parameter choice, namely, for higher values of the 
weakening parameter. The reentrant spinodal scenario 
was one of the conjectures invoked to explain thermody- 
namic anomalies of liquid water, and in particular of the 
divergent-like behavior of response functions in the su- 
percooled regim o 48 ! 49 , even if at the moment the second 
critical point scenario is believed to be more realistic 3 :. 

For the present parameter choice, in our model the 
critical point lies just above the Kauzmann line, while 
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FIG. 5: Response functions at constant pressure (P/rj = 
0.015, 0.030) for pure liquid water as a function of temperature 
(T/rj), for e W w/?7 = 0.3, w — 20, and c w = 0.5. From top to 
bottom, we show the isobaric thermal expansion coefficient 
(r/ap), the isothermal compressibility (t/kt), and the specific 
heat (cp). Numerals beside each plot denote pressure values. 



most of the high-low density liquid transition lies in the 
negative entropy region. This is consistent with the ex- 
perimental observation of two different forms of amorfous 
ice in this region, while the "underlying" liquid phase is 
just an extrapolation of the equation of state for the liq- 
uid. The Kauzmann line displays a cusp (actually a slight 
discontinuity), while intersecting the high-low density liq- 
uid transition. It is noticeable that a similar feature has 
been predicted also by a recent analysis of the poten- 
tial energy landscape of simulated water, performed by 
Sciortino and coworkers^, on the basis of the inherent 
structure theory. 

Let us also report the density response functions and 
the specific heat of the liquid at constant pressure P/rj = 
0.015,0.030, roughly corresponding to 1/5,2/5 of the 
liquid-vapor critical pressure. Also for these calculations, 



details are reported in Ref. 41, We find anomalous be- 
havior, qualitatively similar to that of real liquid water. 
The first response function we consider is the thermal 
expansion coefficient ap = (—dhip/dT)p, which, from 
statistical mechanics, is known to be proportional to the 
entropy-volume cross-correlation. For ordinary fluids, 
ap is always positive, i.e., the local entropy and the local 
specific volume are positively correlated. On the con- 
trary, for our model ap (Fig. [5] top panel) is anomalous. 
As temperature is lowered, the expansion coefficient van- 
ishes (at the temperature of maximum density), and then 
becomes negative. Of course, we do not observe a really 
divergent behavior of this coefficient, but a pronounced 
peak instead. Let us notice that, upon increasing pres- 
sure, the peak is observed to become broader, indicating 
that the liquid is becoming more normal, in agreement 
with experiments. The trend of the isothermal compress- 
ibility kt = (d In p/8P)t is also anomalous (Fig. [5] mid- 
dle panel). For a typical liquid, kt decreases as one 
lowers temperature, because it is proportional to den- 
sity fluctuations, whose magnitude decreases upon de- 
creasing temperature. On the contrary, we observe that 
kt, once reached a minimum, begins to increase upon 
decreasing temperature. The constant pressure specific 
heat cp = (Tds /dT) p (Fig. bottom panel) displays 
a completely analogous behavior, with the minimum oc- 
curring at a higher temperature. 

B. Solution properties 

Let us now consider an ideal inert solute molecule, with 
no Van der Waals interaction with water (e ws = 0), nor 
with other solute molecules (e S5 = 0), and no weaken- 
ing effect on H bonds (c s = 0). As a first analysis of 
the model mixture of water with this kind of (hydropho- 
bic) solute, let us investigate phase diagrams at constant 
temperature and constant pressure, corresponding to the 
"cuts" reported in Fig. 0] We take into account, as a 
composition variable, the solute molar fraction, defined 
as 



where p w and p 5 are determined by Eq. (|14|> . In Fig. 
(top panel), we report a constant temperature phase di- 
agram, computed at T/r) = 0.34. We can observe a first 
order transition between a water-rich and a solute-rich 
phase, which arises continuously from the vapor-liquid 
transition of neat water, as solute concentration is in- 
creased. To determine this transition, we have fixed sev- 
eral different values of water chemical potential /i w , then 
we have determined numerically the value of solute chem- 
ical potential /i s , for which both phases had the same 
pressure P = —id. The phase-separated (coexistence) re- 
gion occupies large part of the diagram, that is, the molar 
fraction of solute which can be dissolved into water, with- 
out giving rise to phase separation, is very small. The 
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FIG. 6: Constant temperature (T/r] = 0.34, top panel) and 
constant pressure (P/rj = 0.015, bottom panel) phase dia- 
grams for a mixture of water (e W w/?7 = 0.3, c w = 0.5, w = 20) 
with an ideal inert solute (e ws = e ss = 0, Cs = 0). Thick lines 
denote pure phase boundaries; thin lines connect coexisting 
phases and denote phase-separated regions. 



behavior is very similar to that of an ordinary solution of 
two disaffine chemical species (at a temperature higher 
than the critical one for the solute), with no peculiar 
anomaly related to the hydrophobic effect. Notice that 
actually we cannot observe any phase transition for neat 
solute (x s — 1), because we have described it as a perfect 
gas, with e ss ~ 0. We have also computed a constant 
pressure phase diagram at P/r) = 0.015, which we have 
reported in Fig. (bottom panel). Here, to compute the 
transition, we have adjusted numerically both chemical 
potentials p w and /i s to impose equality between pres- 
sures of both phases and the reference one. As expected, 
also in this case we observe no phase transition for pure 
solute, whereas a clearly anomalous behavior is observed 
for the water-rich phase boundary. Upon decreasing tem- 
perature, near T/r] « 0.31, the slope of this curve begins 
to change rapidly. Such a behavior gives rise to an ab- 
sorption coefficient x™ /xl (the superscripts denoting the 
water-rich and the solute-rich phase, respectively) with 
a minimum around T/r] w 0.34, well above the temper- 
ature of maximum density at that pressure (see Fig. [SJ. 
This is typical signature of hydrophobic effect. 

Let us now turn to the transfer properties of a sin- 
gle molecule in water, that is, to a dilute solution. Large 
amounts of experimental data are available for this case^. 
According to the Ben-Nairn standard 9 , a transfer process 



(for a chemical species x) can be characterized by means 
of the pseudo-chemical potential p* , related to the ordi- 
nary chemical potential p x by the formula 



p x = p* +T\ogp x . 



(16) 



The use of pseudo-chemical potentials is meant to re- 
move translational entropy contributions, which are not 
directly related to the solvation process. The solvation 
free energy per molecule Ag x , that is, the free energy of 
transfer for a molecule x from the gas phase to the liquid 
phase, can then be defined as 



(17) 



where p* 1 and p* 9 denote pseudo-chemical potentials in 
the liquid and gas phase, respectively. If the gas and liq- 
uid phases coexist in equilibrium, the ordinary chemical 
potentials for the given species must be equal, so that we 
obtain 



a 5 ; 



-Tin 



pi 



(18) 



where p l x and p 9 denote densities for the x species in 
the two phases. Derived quantities, of interest in ex- 
periments, are the solvation entropy 



the solvation enthalpy 



dAg x 



dT 



Ah* x =Ag* + TAs*, 
and the solvation heat-capacity 



OAK 



dT 



(19) 



(20) 



(21) 



In principle, we should distinguish between derivatives 
taken at constant pressure (as stated by definition) or 
along the liquid-vapor equilibrium curve. In particular, 
we could not even use Eq. Ijl8(l . because we would move 
out of the equilibrium curve. Nevertheless, we have veri- 
fied that the difference between the two sets of results is 
negligible, in agreement with experimental observations 9 -, 
and one can usually take the "equilibrium" derivative 
without further care. 

Let us start studying solvation properties for the ideal 
inert solute, in the framework of our model. Water pa- 
rameters are fixed as in the previous case. The temper- 
ature trends of the free energy, entropy, and enthalpy of 
transfer are given in Fig. |7Ji; the transfer heat capacity in 
Fig. 0;. In order to compare with experimental dataS^, 
also reported in Fig. [7J3 and 0i, respectively, all quan- 
tities are evaluated at liquid-vapor coexistence, and for 
very low solute density with respect to water density (di- 
lute solution limit). In practical calculations, we have 
adjusted numerically water and solute chemical poten- 
tials, in order to impose the equilibrium condition (equal 
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FIG. 7: (a) Transfer energy functions (E/rj) vs temperature (T/rj) for an ideal inert molecule in water at liquid-vapor 
coexistence: E — Ag* (solid line), E — TAs* (dashed line), and E = Ah* (dotted line), (b) Corresponding experimental 
data for solvation of argon into watei-. (c) Transfer heat capacity (Acp*) vs temperature (T/rj) for an ideal inert molecule, 
(d) Corresponding experimental data for argon-. Thin lines in (a) and (c) denote transfer quantities for nonbonding water. 



pressure) between liquid and vapor, and to fix the so- 
lute molar fraction. Nevertheless, we have also verified 
that, only for the perfectly inert solute, concentration 
does not affect the results at all, so that we could also set 
an arbitrary value for the solute chemical potential. As 
shown in the previous section, our parameter set for pure 
water corresponds to a liquid-vapor critical temperature 
T/rj w 0.52, and to a temperature of maximum density 
for the liquid phase around T/rj « 0.32, at low pressure. 
Therefore, in order to represent roughly the experimental 
temperature range (between 0° C and 300° C) we report 
model results between T/rj — 0.31 (just below the tem- 
perature of maximum density for pure liquid water) and 
T/rj = 0.40 (about half way between the previous tem- 
perature and the critical temperature). Remarkably, it 
turns out that the model displays the defining features of 
hydrophobic solvation. The solvation free energy is pos- 
itive and large, while the solvation entropy and enthalpy 
are negative at low temperatures and become positive 
upon increasing temperature. The solvation heat capac- 
ity is positive and large, and also the decreasing trend 
with temperature is basically reproduced. A slightly 
increasing trend at high temperature is related to the 
the fact that we are approaching the liquid- vapor critical 
point. Negative solvation entropy at low (room) temper- 
ature is a clear indication that solute insertion into the 
mixture orders the system. The corresponding positive 
(unfavorable) contribution to free energy compensates a 



negative (favorable) enthalpic contribution, giving rise 
to a positive solvation free energy. At higher tempera- 
ture, enthalpic and entropic contributions change sign, 
but they still have the same compensating trend. The 
model also predicts two different temperatures at which 
the transfer enthalpy and entropy vanish (see Fig. [7Jt), 
as observed in experiments (Fig. |?t>). The whole ob- 
served behavior is to be ascribed to the thermodynamics 
of H bonding and, in order to rationalize this fact in the 
model framework, we have also analysed transfer quanti- 
ties, upon removing H bond interactions (see Fig. [7^0:). 
As reasonable, the results are quite similar in the high 
temperature regime, where there is a high probability 
that H bonds are broken by thermal fluctuations, whereas 
they change more and more dramatically upon decreas- 
ing temperature, and, in particular, the regions of nega- 
tive transfer entropy and enthalpy completely disappear. 
This facts confirm that H bonding is the key element for 
system ordering, upon insertion of an inert molecule. Ac- 
cordingly, also the increasing trend of the heat capacity 
upon decreasing temperature is suppressed. The process 
is now dominated by enthalpy, with a large and posi- 
tive transfer free energy (but without a maximum), and 
a positive transfer entropy. Transfer quantities now be- 
have qualitatively as observed in solvation experiments 
of noble gas molecules in ordinary liquids*^, and are 
relatively independent of temperature. Actually, a water 
molecule for which H bond formation has been "turned 
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FIG. 8: (a) Transfer energy functions (E/rj) vs temperature {T/rj) for a water molecule into pure liquid water at liquid-vapor 
coexistence: E = Ag^ (solid line), E — TAs^ (dashed line), and E = Ah^ (dotted line), (b) Corresponding experimental 
data£. Thin lines in (a) denote transfer energies for nonbonding water. 



off" , can be viewed as a nonpolar molecule with only Van 
der Waals interaction energy e ww . 



Let us now consider also the solvation of water in 
its own pure liquid. Corresponding transfer quantities 
obtained by the model are displayed in Fig. |BK, where 
we have reduced the temperature interval, in order to 
compare with available experimental results^, reported 
in Fig. ISJd. With respect to the inert molecule case, 
here absolute values of solvation free energy and entropy 
are considerably smaller. Enthalpy, rather than entropy, 
dominates the solvation process, while all quantities are 
relatively independent of temperature. These features 
characterize a regular transfer process, like the solvation 
of an ordinary fluid molecule from a gas phase into its 
pure liquid phase. In this case, upon removing H bond 
interactions (thin lines in Fig.[HJi), very little changes are 
observed. Let us discuss two issues about these results. 
First, the fact that so little changes are caused by turn- 
ing on or off H bonds can be rationalized on the basis of 
the microscopic model interactions. Insertion of a water 
molecule into pure liquid water should imply in princi- 
ple the formation of new H bonds, but the model is such 
that insertion of a new water molecule also weakens other 
H bonds in its neighborhood, and the two effects nearly 
compensate each other. Second, let us notice that sol- 
vation enthalpy decreases upon increasing temperature, 
that is, the solvation heat capacity is negative, in con- 
trast with experiments. We do not have an explanation 
for this fact, but we can observe that it is basically un- 
changed when H bonds are turned off, that is, when the 
model is reduced to describe a "regular" solvation pro- 
cess. This suggest that there is probably a limitation 
of the lattice description, that anyway has nothing to do 
with peculiarities of water. Indeed, the effect is quantita- 
tively small, so that it is hidden by other large (enthalpic 
and entropic) effects, observed in the case of hydrophobic 
solvation. 



C. "Nonideal" solute and entropy convergence 

So far, we have always turned off all interactions involv- 
ing solute molecules, except excluded volume. Here we 
report some results concerning the role of nonzero solute- 
water interaction parameters (e ws , c s ), still assuming that 
solute molecules have no relevant interaction with one 
another (e ss = 0). In order to have a single parame- 
ter to be varied, we have performed our investigation for 
increasing values of the solute weakening parameter c s , 
and defined solute-water interaction e ws according to the 
following proportionality condition 

At the beginning, we took this assumption as a simple 
trial, but the results we obtained were quite interesting, 
so that we have carried on with the analysis. As far as 
transfer free energies are concerned, we have observed 
a qualitatively unchanged behavior, with a broad maxi- 
mum at some temperature, and free energy values getting 
smaller and smaller, upon increasing e ws and c 5 . On the 
contrary, we have observed peculiar features concerning 
entropies, actually related to one another, which are dis- 
played in Fig. EK- Still upon increasing e W5 and c s , the 
temperature of zero entropy is progressively shifted to- 
wards higher values, while different entropy curves con- 
verge in a very narrow temperature range, relatively close 
to zero entropy temperatures. Moreover, entropy values 
at convergence are negative and relatively small. As one 
can observe in Fig.^, all of these features correspond re- 
markably well to phenomenology observed for the series 
of noble gases^, in particular the entropy convergence, 
which has attracted some interest, due to the fact that a 
similar effect has been observed for the entropies of pro- 
tein unfolding 52 . Because of this unexpected result, we 
have tried and justified a posteriori the working hypothe- 
sis (|22H , and actually a naive explanation may be the fol- 
lowing one. As a first approximation, different hydropho- 
bic species, such as noble gases, may be viewed as hard 
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FIG. 9: (a) Transfer entropy (As*) vs temperature (T/rj) for a solute molecule into pure liquid water at liquid-vapor coexistence, 
for different solute parameter values, with e ss = and e W s/eww = c s /cw. The arrow denotes increasing values of c s and e ws . 
(b) Corresponding experimental data for noble gasea^. 



spheres distinguished by their diameter, that is their vol- 
ume, only. In this way, the proportionality condition can 
be conceived as a trick that, in the framework of the lat- 
tice model, mimics the fact that "a fraction" c s /c w of the 
site occupied by the solute actually "behaves like water" . 
As a consequence, higher values of the ratio c s /c w should 
correspond to smaller solute molecules, which turns out 
to be consistent, comparing Figs.|HK and|5jD. Let us no- 
tice also that, for parameter values as in Fig. the 
"fractions that behaves like solute" (1 — c s /c w ) correlate 
well with the squares of the hard sphere diameters of the 
corresponding substances^, which is consistent with the 
common assumption that hydrophobicity is proportional 
to exposed surface. 

V. DISCUSSION AND CONCLUSIONS 

In this paper we have considered a 3-dimensional lat- 
tice fluid model of water, which we had previously shown 
to exhibit realistic thermodynamic anomalies^, and ex- 
tended the model to describe aqueous solutions. The 
motivation for this work resides mainly in the interest 
for the hydrophobic effect, whose relevance for biological 
processes such as protein folding, taking place in aqueous 
solutions, has been more and more recognized in the lat- 
est years. Moreover, a simplified but accurate modelling 
of water is an appealing issue, in view of investigations 
on such processes, because detailed water models may be 
extremely time consuming^. 

As far as water is concerned, our model is a simplified 
version of a previous model proposed by Roberts and 
Debenedetti^, without a distinction between hydrogen 
bond donors and acceptors. In the framework of this 
model, the microscopic description of water anomalies, is 
essentially based on the competition between an isotropic 
(Van der Waals like) interaction and an highly directional 
(H bonding) interaction, and on the difference between 
the respective optimal interaction distances. In the lat- 



tice environment, the latter is taken into account by a 
trick, that is, the weakening effect of a water molecule 
near a formed bond. The same assumptions also ac- 
counts for possible perturbations of the electronic den- 
sity, due to interaction with other water molecules. We 
have extended the weakening trick to take into account 
the presence of a different chemical species (solute) , with 
no internal degrees of freedom (bonding arms). Calcu- 
lations have been performed in a generalized first-order 
approximation on a tetrahedral cluster, which requires 
small computational effort, and had been shown to be 
quite accurate for the pure water model^i. 

Due to the fact that we have chosen slightly differ- 
ent interaction parameters for the present investigation, 
we have first analysed again pure water behavior. At 
constant pressure, the typical thermodynamic anomalies 
are reproduced, with a density maximum, and a mini- 
mum of isothermal compressibility and specific heat. In 
the ordinary temperature and pressure region, the tem- 
perature of maximum density decreases upon increasing 
pressure, as observed in experiments. Moreover, as far 
as the supercooled regime is concerned, there is still evi- 
dence of a second (metastable) critical point, which ter- 
minates a line of coexistence between two liquid phases 
at different densities. Let us recall that the pure water 
model could predict, for different values of the weaken- 
ing parameter c w , two different scenarios, that is, with 
or without a reentrant spinodal. The present parameter 
choice predicts a nonreentrant spinodal. The reentrant 
spinodal scenario was the first conjecture put forth to 
justify water anomalies, whereas the most recent and ac- 
curate molecular dynamics simulations of water suggest 
a scenario with a nonreentrant spinodal and a metastable 
liquid-liquid critical pointi As far as the critical point 
is concerned, in our calculation, it lies at some temper- 
ature just above the Kauzmann line, at which the con- 
figurational entropy vanishes, while the Kauzmann line 
displays a cusp upon crossing the metastable coexistence 
line. All of these features turn out to be in a remarkably 
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good agreement with the inherent state analysis of the 
potential energy landscape of simulated water, recently 
performed by Sciortino and coworkers^. 

As far as the solution model is concerned, we have 
mainly considered solutions of inert, that is, hydrophobic 
solutes. First of all, we have investigated phase diagrams 
for arbitrary solute concentration, pointing out a mini- 
mum of solubility as a function of temperature, as typical 
for hydrophobic solutions. We have investigated in more 
detail the dilute solution limit, at liquid- vapor equilib- 
rium, for which many experimental data are available. 
Solvation quantities turn out to exhibit peculiar features 
that are believed to be the fingerprints of hydrophobicity. 
The solvation free energy is positive (unfavorable solva- 
tion), while entropy and enthalpy are negative at low 
temperatures and positive at high temperatures. The 
solvation heat capacity is large and decreases upon in- 
creasing temperature. These results compare qualita- 
tively well with solvation experiments for noble gases in 
water. Let us notice that a previous lattice model by 
Besseling and Lyklema™, based on a different descrip- 
tion of water interactions, was also able to account for 
the qualitative behavior of free energy, enthalphy, and 
entropy of transfer. Nevertheless, it failed in reproduc- 
ing the correct temperature trend of the transfer heat 
capacity, which, according to some authors^ is a key fea- 
ture, revealing the cooperative nature of the hydrophobic 
effect. Let us recall, by the way, that the same difficulty 
about heat capacity is encountered by the information 
theory approach^*. 

We have investigated explicitly on the effect of H bond- 
ing, in the framework of our model, performing calcu- 
lations also when this interaction is completely turned 
off. In this case, we have obtained transfer quantities 
that approach the ones computed with H bonds at high 
temperatures, but that largely deviates from them upon 
decreasing temperature, that is, in the region were H 
bonding begins to dominate. In particular, we have ob- 
served that, while disaffinity between solute and solvent 
is left (the solvation free energy is still positive), this is 
mainly of enthalpic nature. Both the enthalpy and en- 
tropy of solvation remain positive at all temperatures, so 
that also the typical strong temperature dependence of 
hydrophobic solvation, disappears. 

We have also taken into account the solvation process 
of water into its own pure liquid, for which experimental 
data are available. We have found qualitative agreement, 
as far as the values of transfer free energy, entropy and en- 
thalpy are concerned, but we have observed some discrep- 
ancy in the temperature dependence of enthalpy, indicat- 
ing a negative solvation heat capacity, in disagreement 
with experiments. We have verified that the same kind 
of discrepancy can be observed if H bonding is turned off, 
that is, for an ordinary lattice gas. Therefore, we sug- 
gest that the discrepancy is to be related to an intrinsic 
limitation of the lattice environment, that has nothing 
to do with peculiarities of the water model. The effect is 



relatively small, so that it is completely invisible, when 
the dominant effect of H bonds is introduced. 

Let us notice that the results concerning transfer quan- 
tities are qualitatively similar to those observed for a 2- 
dimensional lattice model recently investigated by us^, 
which in turn is a simplified "lattice version" of the 
Mercedes Benz model investigated by Dill and cowork- 
ers. In spite of similarities in the experimental tempera- 
ture range, there is one major difference between the two 
models. The 2-dimensional model is not able to account 
for the metastable critical point of water, probably due 
to impossibility of reproducing a high density ordered 
packing of water molecules. As a consequence, thermo- 
dynamic anomalies are entirely due to the presence of a 
reentrant spinodal, which is no longer believed to be a re- 
alistic scenario for real water—. In this sense, the present 
3-dimensional model seems to provide a more coherent 
view of water thermodynamics, relating each other neat 
water anomalies and hydrophobic effect. 

Moreover, we have shown that, with quite a reasonable 
assumption for the solute interaction parameters e ws ,c s 
(attempting to describe solutes of different volume in the 
lattice framework), the model is able to reproduce also 
the entropy convergence phenomenon, in a qualitatively 
correct way. Such a phenomenon has been theoretically 
described, for instance, by the simplified molecular the- 
ory of Debenedetti and coworkers^. Moreover, an almost 
quantitative explanation has been proposed by Pratt and 
coworkers^, on the basis of the information theory ap- 
proach, which we mentioned in the Introduction. In the 
cited work, the authors argue that entropy convergence, 
and in particular the negative entropy at convergence, 
are related to the weak temperature dependence of free 
volume fluctuations in liquid water, that is, isothermal 
compressibility. Although we cannot provide a clear ex- 
planation of why our lattice model, with the proportion- 
ality assumption l)22|). exhibits qualitatively correct be- 
havior, let us notice that such result is somehow consis- 
tent with the previous explanation. In fact, entropy con- 
vergence occurs very close to the minimum of isothermal 
compressibility, as one can argue from Fig. [5] that is, in 
a region where compressibility is nearly constant, as a 
function of temperature. Nonetheless, the proportional- 
ity assumption IL'L'I) still remains not well justified. 

Let us finally recall that our model, at least in the 
present treatment, is not able to provide microscopic 
structural details as simulations do, but its most appeal- 
ing feature is simplicity. In the present paper we have 
shown that, in spite of this, the model yields a qualita- 
tively coherent description of peculiar thermodynamics of 
water, not only as a pure substance but also as a solvent, 
and is consistent with predictions based on much more 
sophisticated models and simulations. Moreover, thanks 
to the 3-dimensional embedding, it may be suitable for 
quite a realistic analysis of more complex, for instance 
polymeric, solutes. We are going to report about such 
investigations in a forthcoming article. 
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